{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Graphical Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets import fetch_openml\n",
    "from prml import bayesnet as bn\n",
    "\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = bn.discrete([0.1, 0.9])\n",
    "f = bn.discrete([0.1, 0.9])\n",
    "\n",
    "g = bn.discrete([[[0.9, 0.8], [0.8, 0.2]], [[0.1, 0.2], [0.2, 0.8]]], b, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(proba=[0.1 0.9])\n",
      "f: DiscreteVariable(proba=[0.1 0.9])\n",
      "g: DiscreteVariable(proba=[0.315 0.685])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "g.observe(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
      "f: DiscreteVariable(proba=[0.25714286 0.74285714])\n",
      "g: DiscreteVariable(observed=[1. 0.])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.observe(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b: DiscreteVariable(observed=[1. 0.])\n",
      "f: DiscreteVariable(proba=[0.11111111 0.88888889])\n",
      "g: DiscreteVariable(observed=[1. 0.])\n"
     ]
    }
   ],
   "source": [
    "print(\"b:\", b)\n",
    "print(\"f:\", f)\n",
    "print(\"g:\", g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.3 Illustration: Image de-noising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fa5a28f1a90>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAALNElEQVR4nO3dT6hm9X3H8fenJtkYoWOlwzAxNS3usjBFXEmxiwTrZsxG4mpCCjeLWtJdJFlECIFQ2nRZMEQyLakhoNZBShMrIWYVHMXqqCTaMJIZxhlkWmpWafTbxT0jN+P9N895znOeO9/3Cx6e5zn3ued8Pd7P/H7n97vn/lJVSLr2/d7cBUhaDcMuNWHYpSYMu9SEYZea+NAqD5bEoX9pYlWV7baPatmT3J3k50neSPLgmH1JmlYWnWdPch3wC+DTwFngOeD+qnp1l++xZZcmNkXLfgfwRlX9sqp+A3wfODZif5ImNCbsR4FfbXl/dtj2O5JsJDmV5NSIY0kaafIBuqp6GHgY7MZLcxrTsp8Dbt7y/mPDNklraEzYnwNuTfKJJB8BPgecXE5ZkpZt4W58Vf02yQPAD4HrgEeq6pWlVSZpqRaeelvoYF6zS5Ob5JdqJB0chl1qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhOGXWrCsEtNGHapCcMuNWHYpSYMu9SEYZeaMOxSE4ZdasKwS00YdqkJwy41YdilJgy71IRhl5ow7FIThl1qwrBLTRh2qQnDLjWx8JLN0lhjVxBOtl2sdGn7H3PsdTQq7EnOAO8A7wK/rarbl1GUpOVbRsv+51X19hL2I2lCXrNLTYwNewE/SvJ8ko3tPpBkI8mpJKdGHkvSCBkziJHkaFWdS/KHwNPAX1fVs7t8froREx04DtBNo6q2LW5Uy15V54bni8ATwB1j9idpOguHPcn1SW64/Br4DHB6WYVJWq4xo/GHgSeG7syHgH+pqn9fSlW6KlN2V9dZ1//uRY26Zr/qg3nNPgl/6Fev3TW7pIPDsEtNGHapCcMuNWHYpSa8xXUFHC1fzDqPeB9EtuxSE4ZdasKwS00YdqkJwy41YdilJgy71IRhl5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmDLvUhPezr8DYlUuu1ZVPvM9/tWzZpSYMu9SEYZeaMOxSE4ZdasKwS00YdqkJ59nXwJx/H73rsTvas2VP8kiSi0lOb9l2Y5Knk7w+PB+atkxJY+2nG/9d4O4rtj0IPFNVtwLPDO8lrbE9w15VzwKXrth8DDgxvD4B3LvkuiQt2aLX7Ier6vzw+i3g8E4fTLIBbCx4HElLMnqArqoqyY53NFTVw8DDALt9TtK0Fp16u5DkCMDwfHF5JUmawqJhPwkcH14fB55cTjmSppJ93Ev9KHAXcBNwAfga8K/AD4CPA28C91XVlYN42+3LbvwEDur97JpGVW37P23PsC+TYZ+GYddWO4XdX5eVmjDsUhOGXWrCsEtNGHapCW9xvQbsNmLun2vWZbbsUhOGXWrCsEtNGHapCcMuNWHYpSYMu9SE8+zXuLHLPY+dp/euufVhyy41YdilJgy71IRhl5ow7FIThl1qwrBLTTjP3tzYefi97Pb9zsGvli271IRhl5ow7FIThl1qwrBLTRh2qQnDLjXhPLt2NeU8vPfKr9aeLXuSR5JcTHJ6y7aHkpxL8uLwuGfaMiWNtZ9u/HeBu7fZ/g9Vddvw+LflliVp2fYMe1U9C1xaQS2SJjRmgO6BJC8N3fxDO30oyUaSU0lOjTiWpJGyn0GSJLcAT1XVJ4f3h4G3gQK+Dhypqi/sYz+uMniNmXPhSAfotldV256YhVr2qrpQVe9W1XvAt4E7xhQnaXoLhT3JkS1vPwuc3umzktbDnvPsSR4F7gJuSnIW+BpwV5Lb2OzGnwG+OGGNWmNjutJT3isPdvOvtK9r9qUdzGt2bTH1z17XsC/1ml3SwWPYpSYMu9SEYZeaMOxSE4ZdasKwS00YdqkJwy41YdilJgy71IRhl5ow7FIT/ilpjTLnX6rR1bFll5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmnGdvznnyPmzZpSYMu9SEYZeaMOxSE4ZdasKwS00YdqkJ59mvcQd5Hr3rKqxT2bNlT3Jzkh8neTXJK0m+NGy/McnTSV4fng9NX66kRe25PnuSI8CRqnohyQ3A88C9wOeBS1X1zSQPAoeq6st77OvgNjMHlC17Pwuvz15V56vqheH1O8BrwFHgGHBi+NgJNv8BkLSmruqaPcktwKeAnwGHq+r88KW3gMM7fM8GsLF4iZKWYc9u/PsfTD4K/AT4RlU9nuR/qur3t3z9v6tq1+t2u/GrZze+n4W78QBJPgw8Bnyvqh4fNl8YrucvX9dfXEahkqaxn9H4AN8BXquqb2350kng+PD6OPDk8ssTbLbOiz7mlmThh5ZrP6PxdwI/BV4G3hs2f4XN6/YfAB8H3gTuq6pLe+xr/p++A2gdQrsoQ7t6O3Xj933NvgyGfTGGXVdj1DW7pIPPsEtNGHapCcMuNWHYpSa8xXUJDvJo+V4cTb922LJLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhPOsw+u5bny3TiP3octu9SEYZeaMOxSE4ZdasKwS00YdqkJwy410Wae/VqeR3euXPthyy41YdilJgy71IRhl5ow7FIThl1qwrBLTexnffabk/w4yatJXknypWH7Q0nOJXlxeNwzfbmLG7NO+Lo/pP3Yz/rsR4AjVfVCkhuA54F7gfuAX1fV3+37YC7ZLE1upyWb9/wNuqo6D5wfXr+T5DXg6HLLkzS1q7pmT3IL8CngZ8OmB5K8lOSRJId2+J6NJKeSnBpVqaRR9uzGv//B5KPAT4BvVNXjSQ4DbwMFfJ3Nrv4X9tiH3XhpYjt14/cV9iQfBp4CflhV39rm67cAT1XVJ/fYj2GXJrZT2PczGh/gO8BrW4M+DNxd9lng9NgiJU1nP6PxdwI/BV4G3hs2fwW4H7iNzW78GeCLw2DebvuyZZcmNqobvyyGXZrewt14SdcGwy41YdilJgy71IRhl5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhOrXrL5beDNLe9vGrato3WtbV3rAmtb1DJr+6OdvrDS+9k/cPDkVFXdPlsBu1jX2ta1LrC2Ra2qNrvxUhOGXWpi7rA/PPPxd7Outa1rXWBti1pJbbNes0tanblbdkkrYtilJmYJe5K7k/w8yRtJHpyjhp0kOZPk5WEZ6lnXpxvW0LuY5PSWbTcmeTrJ68PztmvszVTbWizjvcsy47Oeu7mXP1/5NXuS64BfAJ8GzgLPAfdX1asrLWQHSc4At1fV7L+AkeTPgF8D/3R5aa0kfwtcqqpvDv9QHqqqL69JbQ9xlct4T1TbTsuMf54Zz90ylz9fxBwt+x3AG1X1y6r6DfB94NgMday9qnoWuHTF5mPAieH1CTZ/WFZuh9rWQlWdr6oXhtfvAJeXGZ/13O1S10rMEfajwK+2vD/Leq33XsCPkjyfZGPuYrZxeMsyW28Bh+csZht7LuO9SlcsM742526R5c/HcoDug+6sqj8F/gL4q6G7upZq8xpsneZO/xH4EzbXADwP/P2cxQzLjD8G/E1V/e/Wr8157rapayXnbY6wnwNu3vL+Y8O2tVBV54bni8ATbF52rJMLl1fQHZ4vzlzP+6rqQlW9W1XvAd9mxnM3LDP+GPC9qnp82Dz7uduurlWdtznC/hxwa5JPJPkI8Dng5Ax1fECS64eBE5JcD3yG9VuK+iRwfHh9HHhyxlp+x7os473TMuPMfO5mX/68qlb+AO5hc0T+v4CvzlHDDnX9MfCfw+OVuWsDHmWzW/d/bI5t/CXwB8AzwOvAfwA3rlFt/8zm0t4vsRmsIzPVdiebXfSXgBeHxz1zn7td6lrJefPXZaUmHKCTmjDsUhOGXWrCsEtNGHapCcMuNWHYpSb+H6RpBIl+5K8zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mnist = fetch_openml(\"mnist_784\")\n",
    "x = mnist.data[0]\n",
    "binarized_img = (x > 127).astype(np.int).reshape(28, 28)\n",
    "plt.imshow(binarized_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fa598be71d0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAMsklEQVR4nO3dT6hc5R3G8eep2o0KTWp7ucTY2OLOhfZesgrFLpQ0m+hGdBWx9Lqoxe4UuzAgQiitxVUh1mAsVhGMNYhUUxHjSnIT0pg/1FiJmBBzlbQ0rqzm18WcyDXOn5s558w5Z37fDwwz98zcOb97Zp573vO+M+d1RAjA9PtW0wUAmAzCDiRB2IEkCDuQBGEHkrh8kiuzPZVd/3Nzc0Pv379/f2Prr3vdXdX0a1aniHC/5S4z9GZ7o6QnJF0m6U8RsW3E46cy7KO2od13209k/XWvu6uafs3qVHnYbV8m6T1Jt0o6KWmfpLsj4uiQ3yHsE15/l9+0dWr6NavToLCXOWZfL+n9iPggIj6X9LykzSWeD0CNyoR9jaSPlv18slj2NbYXbC/aXiyxLgAl1d5BFxHbJW2XprcZD3RBmT37KUlrl/18bbEMQAuVCfs+STfYvt72tyXdJWl3NWUBqNrYzfiI+ML2/ZJeU2/obUdEHKmssg5puud22Pqnude5jIx/d6lx9kteGcfsE0fY86lj6A1AhxB2IAnCDiRB2IEkCDuQBGEHkpjo99mbVHaIsc4hqjqHxxhamz7D3i/z8/MD72PPDiRB2IEkCDuQBGEHkiDsQBKEHUgizdDbqCGoJie4ZHgMl2Lc9wt7diAJwg4kQdiBJAg7kARhB5Ig7EAShB1IIs04+yjTOtbd9Fd765x0ssxXgzOedZc9O5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kMTXj7HWPm5YZL25yTLfL48VlPyNQ5vencRy+VNhtn5B0TtKXkr6IiMEnrQbQqCr27D+NiE8reB4ANeKYHUiibNhD0uu299te6PcA2wu2F20vllwXgBJcphPD9pqIOGX7+5L2SPpVROwd8vjazupIB1071bndmtTm1ywi+hZXas8eEaeK6yVJL0laX+b5ANRn7LDbvtL21RduS7pN0uGqCgNQrTK98TOSXiqaM5dL+ktE/K1MMW1u7ra82dZ0CWPpat1dVeqY/ZJXNuKYfVpPRtBkfwLq0fJ//tUfswPoDsIOJEHYgSQIO5AEYQeSaNVXXMv0cLa5dxTj4TWtFnt2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiiVePsZbT5W2+jdLl2dAd7diAJwg4kQdiBJAg7kARhB5Ig7EAShB1IYmrG2duMmU/6a/PfPY3YswNJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEp0aZx82Ltvm73y3+fvqda+7q69Z3ZrYLiP37LZ32F6yfXjZstW299g+XlyvqqU6AJVZSTP+aUkbL1r2kKQ3IuIGSW8UPwNosZFhj4i9ks5etHizpJ3F7Z2Sbq+4LgAVG/eYfSYiThe3P5Y0M+iBthckLYy5HgAVKd1BFxFhe2BvQ0Rsl7RdkoY9DkC9xh16O2N7VpKK66XqSgJQh3HDvlvSluL2FkkvV1MOgLp4BWPAz0m6RdI1ks5IekTSXyW9IOk6SR9KujMiLu7E6/dcNONbps2fARilzPfh2/x3lRURff+4kWGvEmFvH8I+fQaFnY/LAkkQdiAJwg4kQdiBJAg7kARfcZ2AJnu8mz5dc52vWZtf8zZizw4kQdiBJAg7kARhB5Ig7EAShB1IgrADSXRqnL2ryk7ZXGY8ue7pouusvcvfyKvTsO0yPz8/8D727EAShB1IgrADSRB2IAnCDiRB2IEkCDuQxETPLjs/Px+Li4uDi2FMtnOaPMNrm1/zhs9hwNllgcwIO5AEYQeSIOxAEoQdSIKwA0kQdiAJZnGtQJvHe+vW5Hnpp3m7ljH2OLvtHbaXbB9etmyr7VO2DxaXTVUWC6B6K2nGPy1pY5/lf4iIm4rLq9WWBaBqI8MeEXslnZ1ALQBqVKaD7n7bh4pm/qpBD7K9YHvR9uAPxQOo3Yo66Gyvk/RKRNxY/Dwj6VNJIelRSbMRce8KnocOuilDB137VPpFmIg4ExFfRsR5SU9KWl+mOAD1GyvstmeX/XiHpMODHgugHUaeN972c5JukXSN7ZOSHpF0i+2b1GvGn5B0X401tl7m5mSd5yBAtfhQDRpT9r2X+Z/sMJy8AkiOsANJEHYgCcIOJEHYgSRaNWVz5k+i1WWat+k0/211YM8OJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0m0apy9q+OibR7vrXvdfE11PMO2W12vGXt2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUhiomGfm5tTRAy8jFLX75b9fdtDL11WdrtNq7LbpYn3C3t2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCWVynXJfHwpv8jEKbz1EwytizuNpea/tN20dtH7H9QLF8te09to8X16uqLhpAdUbu2W3PSpqNiAO2r5a0X9Ltku6RdDYittl+SNKqiHhwxHN1dzfTUezZx5Nyzx4RpyPiQHH7nKRjktZI2ixpZ/Gwner9AwDQUpd0Djrb6yTdLOkdSTMRcbq462NJMwN+Z0HSwvglAqjCijvobF8l6S1Jj0XELtv/iYjvLLv/3xEx9LidZvzk0YwfT8pmvCTZvkLSi5KejYhdxeIzxfH8heP6pSoKBVCPkc149/6FPSXpWEQ8vuyu3ZK2SNpWXL9cS4XLNHH63bavu+3avAccpqt1D7OS3vgNkt6W9K6k88Xih9U7bn9B0nWSPpR0Z0ScHfFcpd61bQ0cYR9sGkPTdoOa8Z36UE1bA0fYByPsk1fqmB1A9xF2IAnCDiRB2IEkCDuQRKumbB6lrT27Xe4tH6Wt2xyXjj07kARhB5Ig7EAShB1IgrADSRB2IAnCDiTRqXH2Ok3zWPkwWcfRu3wmmnGxZweSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJNKMs0/zOPqwMeFp/rvLaPM4el2fAWDPDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJrGR+9rWSnpE0IykkbY+IJ2xvlfQLSZ8UD304Il6tq9CyRo1NTuv3m7tad2Z1vWYrmZ99VtJsRBywfbWk/ZJul3SnpM8i4ncrXlnJKZvrNK1hRz6DpmweuWePiNOSThe3z9k+JmlNteUBqNslHbPbXifpZknvFIvut33I9g7bqwb8zoLtRduLpSoFUMrIZvxXD7SvkvSWpMciYpftGUmfqncc/6h6Tf17RzwHzXigZoOa8SsKu+0rJL0i6bWIeLzP/eskvRIRN454HsIO1GxQ2Ec24917lz8l6djyoBcddxfcIelw2SIB1GclvfEbJL0t6V1J54vFD0u6W9JN6jXjT0i6r+jMG/Zcnd2zD1N2r0+ron26/JqUasZXhbCPt+42v7GmVZdfk7Gb8QCmA2EHkiDsQBKEHUiCsANJEHYgiYmeSnpubk6Li/V8RL7sUEiTQyl1Dt0xLDieafy72LMDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKT/orrJ5I+XLboGvVObdVGba2trXVJ1DauKmv7QUR8r98dEw37N1ZuL0bEfGMFDNHW2tpal0Rt45pUbTTjgSQIO5BE02Hf3vD6h2lrbW2tS6K2cU2ktkaP2QFMTtN7dgATQtiBJBoJu+2Ntv9p+33bDzVRwyC2T9h+1/bBpuenK+bQW7J9eNmy1bb32D5eXPedY6+h2rbaPlVsu4O2NzVU21rbb9o+avuI7QeK5Y1uuyF1TWS7TfyY3fZlkt6TdKukk5L2Sbo7Io5OtJABbJ+QNB8RjX8Aw/ZPJH0m6ZkLU2vZ/q2ksxGxrfhHuSoiHmxJbVt1idN411TboGnG71GD267K6c/H0cSefb2k9yPig4j4XNLzkjY3UEfrRcReSWcvWrxZ0s7i9k713iwTN6C2VoiI0xFxoLh9TtKFacYb3XZD6pqIJsK+RtJHy34+qXbN9x6SXre93/ZC08X0MbNsmq2PJc00WUwfI6fxnqSLphlvzbYbZ/rzsuig+6YNEfFjST+T9MuiudpK0TsGa9PY6R8l/Ui9OQBPS/p9k8UU04y/KOnXEfHf5fc1ue361DWR7dZE2E9JWrvs52uLZa0QEaeK6yVJL6l32NEmZy7MoFtcLzVcz1ci4kxEfBkR5yU9qQa3XTHN+IuSno2IXcXixrddv7omtd2aCPs+STfYvt72tyXdJWl3A3V8g+0ri44T2b5S0m1q31TUuyVtKW5vkfRyg7V8TVum8R40zbga3naNT38eERO/SNqkXo/8vyT9pokaBtT1Q0n/KC5Hmq5N0nPqNev+p17fxs8lfVfSG5KOS/q7pNUtqu3P6k3tfUi9YM02VNsG9ZrohyQdLC6bmt52Q+qayHbj47JAEnTQAUkQdiAJwg4kQdiBJAg7kARhB5Ig7EAS/wcn8WtNY83dTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "indices = np.random.choice(binarized_img.size, size=int(binarized_img.size * 0.1), replace=False)\n",
    "noisy_img = np.copy(binarized_img)\n",
    "noisy_img.ravel()[indices] = 1 - noisy_img.ravel()[indices]\n",
    "plt.imshow(noisy_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "markov_random_field = np.array([\n",
    "        [[bn.discrete([0.5, 0.5], name=f\"p(z_({i},{j}))\") for j in range(28)] for i in range(28)], \n",
    "        [[bn.DiscreteVariable(2) for _ in range(28)] for _ in range(28)]])\n",
    "a = 0.9\n",
    "b = 0.9\n",
    "pa = [[a, 1 - a], [1 - a, a]]\n",
    "pb = [[b, 1 - b], [1 - b, b]]\n",
    "for i, j in itertools.product(range(28), range(28)):\n",
    "    bn.discrete(pb, markov_random_field[0, i, j], out=markov_random_field[1, i, j], name=f\"p(x_({i},{j})|z_({i},{j}))\")\n",
    "    if i != 27:\n",
    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i + 1, j]], name=f\"p(z_({i},{j}), z_({i+1},{j}))\")\n",
    "    if j != 27:\n",
    "        bn.discrete(pa, out=[markov_random_field[0, i, j], markov_random_field[0, i, j + 1]], name=f\"p(z_({i},{j}), z_({i},{j+1}))\")\n",
    "    markov_random_field[1, i, j].observe(noisy_img[i, j], proprange=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fa59879cd30>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAALWklEQVR4nO3dT6yldX3H8fenqBsk6VDayQSx2IadCyyEFWnoQkPZgBsiqzE2uS5KY3cSXUhiTEzT2mUTjMRpYzEmQCGkqVJixJVhIBQGiELNEGcyzJRMjbiywreL+wy5wv035znnPM+93/crOTnnPOfcc77zzP3c3+/5/c5zfqkqJB1+vzd1AZLWw7BLTRh2qQnDLjVh2KUmPrDON0vi0L+0YlWV7baPatmT3J7kp0leS3LfmNeStFpZdJ49yRXAz4BPAmeAZ4B7qurlXX7Gll1asVW07LcAr1XVz6vqN8B3gTtHvJ6kFRoT9muBX2y5f2bY9juSbCQ5meTkiPeSNNLKB+iq6gHgAbAbL01pTMt+Frhuy/2PDNskzdCYsD8D3JDkY0k+BHwGeHw5ZUlatoW78VX12yT3At8HrgAerKqXllaZpKVaeOptoTfzmF1auZV8qEbSwWHYpSYMu9SEYZeaMOxSE4ZdamKt57PP2V5TkMm2sxnSgWHLLjVh2KUmDLvUhGGXmjDsUhOGXWrCqbeBU2s67GzZpSYMu9SEYZeaMOxSE4ZdasKwS00YdqkJ59l1aO122nLHz1XYsktNGHapCcMuNWHYpSYMu9SEYZeaMOxSE86za6XWuUrw5Rhb10Gcpx8V9iSngbeAt4HfVtXNyyhK0vIto2X/i6p6cwmvI2mFPGaXmhgb9gJ+kOTZJBvbPSHJRpKTSU6OfC9JI2TMQEWSa6vqbJI/Ap4E/qaqnt7l+fMcrdHKzHWAbqw5D9BV1bbFjWrZq+rscH0BeBS4ZczrSVqdhcOe5MokV126DXwKOLWswiQt15jR+KPAo0N35gPAv1bVfyylqkNm1ctBH9auspZr1DH7Zb9Z02N2w374tDtml3RwGHapCcMuNWHYpSYMu9SEp7hqtuY84n0Q2bJLTRh2qQnDLjVh2KUmDLvUhGGXmjDsUhPOs8/Aqs+Kk8CWXWrDsEtNGHapCcMuNWHYpSYMu9SEYZeacJ59BuY8j77Kb76d87/7MLJll5ow7FIThl1qwrBLTRh2qQnDLjVh2KUmnGdfg87zyZ3/7XOzZ8ue5MEkF5Kc2rLt6iRPJnl1uD6y2jIljbWfbvy3gdvfs+0+4KmqugF4argvacb2DHtVPQ1cfM/mO4ETw+0TwF1LrkvSki16zH60qs4Nt98Aju70xCQbwMaC7yNpSUYP0FVVJdnxbIeqegB4AGC350larUWn3s4nOQYwXF9YXkmSVmHRsD8OHB9uHwceW045klYl+/jO8oeA24BrgPPAV4B/A74HfBR4Hbi7qt47iLfda9mNXzO/k76fqtr2P3XPsC+TYV8/w97PTmH347JSE4ZdasKwS00YdqkJwy414Smuh9xeo+2O1vdhyy41YdilJgy71IRhl5ow7FIThl1qwrBLTTjP3tyq59Fdsnk+bNmlJgy71IRhl5ow7FIThl1qwrBLTRh2qQnn2TXKKr+d2HPtl8uWXWrCsEtNGHapCcMuNWHYpSYMu9SEYZeacJ5dk1nnCsLaR8ue5MEkF5Kc2rLt/iRnkzw/XO5YbZmSxtpPN/7bwO3bbP/HqrpxuPz7csuStGx7hr2qngYurqEWSSs0ZoDu3iQvDN38Izs9KclGkpNJTo54L0kjZT+DJEmuB56oqo8P948CbwIFfBU4VlWf28frOCJzyEw5yOaJMNurqm13zEIte1Wdr6q3q+od4JvALWOKk7R6C4U9ybEtdz8NnNrpuZLmYc959iQPAbcB1yQ5A3wFuC3JjWx2408Dn19hjZqxMV1p59nXa1/H7Et7M4/ZtcXY3z2P2be31GN2SQePYZeaMOxSE4ZdasKwS014iqsOLL9q+vLYsktNGHapCcMuNWHYpSYMu9SEYZeaMOxSE86zaxRPUz04bNmlJgy71IRhl5ow7FIThl1qwrBLTRh2qQnn2bUr59EPD1t2qQnDLjVh2KUmDLvUhGGXmjDsUhOGXWrCefZD7iDPk/u978u1Z8ue5LokP0zycpKXknxh2H51kieTvDpcH1l9uZIWtef67EmOAceq6rkkVwHPAncBnwUuVtXXk9wHHKmqL+7xWge3mTmgbNn7WXh99qo6V1XPDbffAl4BrgXuBE4MTzvB5h8ASTN1WcfsSa4HPgH8BDhaVeeGh94Aju7wMxvAxuIlSlqGPbvx7z4x+TDwI+BrVfVIkl9W1e9vefx/q2rX43a78etnN76fhbvxAEk+CDwMfKeqHhk2nx+O5y8d119YRqGSVmM/o/EBvgW8UlXf2PLQ48Dx4fZx4LHllyfYbJ0XvUwtycIXLdd+RuNvBX4MvAi8M2z+EpvH7d8DPgq8DtxdVRf3eK3pf/sOoDmEdlGGdv126sbv+5h9GQz7Ygy7LseoY3ZJB59hl5ow7FIThl1qwrBLTXiK6z4d5BHxMRxNPzxs2aUmDLvUhGGXmjDsUhOGXWrCsEtNGHapCefZB86j67CzZZeaMOxSE4ZdasKwS00YdqkJwy41YdilJpxnP+ScR9cltuxSE4ZdasKwS00YdqkJwy41YdilJgy71MR+1me/LskPk7yc5KUkXxi235/kbJLnh8sdqy93nsasQb7qi3TJftZnPwYcq6rnklwFPAvcBdwN/Lqq/n7fbzbjJZvHfHmFodKc7LRk856foKuqc8C54fZbSV4Brl1ueZJW7bKO2ZNcD3wC+Mmw6d4kLyR5MMmRHX5mI8nJJCdHVSpplD278e8+Mfkw8CPga1X1SJKjwJtAAV9ls6v/uT1ew268tGI7deP3FfYkHwSeAL5fVd/Y5vHrgSeq6uN7vI5hl1Zsp7DvZzQ+wLeAV7YGfRi4u+TTwKmxRUpanf2Mxt8K/Bh4EXhn2Pwl4B7gRja78aeBzw+Debu91oFt2W29dVCM6sYvi2GXVm/hbrykw8GwS00YdqkJwy41YdilJgy71MRaw37TTTdRVSu5jOWpojrsbNmlJgy71IRhl5ow7FIThl1qwrBLTRh2qYl1n+L6P8DrWzZdw+ZXW83RXGuba11gbYtaZm1/XFV/uN0Daw37+948OVlVN09WwC7mWttc6wJrW9S6arMbLzVh2KUmpg77AxO//27mWttc6wJrW9Raapv0mF3S+kzdsktaE8MuNTFJ2JPcnuSnSV5Lct8UNewkyekkLw7LUE+6Pt2wht6FJKe2bLs6yZNJXh2ut11jb6LaZrGM9y7LjE+676Ze/nztx+xJrgB+BnwSOAM8A9xTVS+vtZAdJDkN3FxVk38AI8mfA78G/vnS0lpJ/g64WFVfH/5QHqmqL86ktvu5zGW8V1TbTsuMf5YJ990ylz9fxBQt+y3Aa1X186r6DfBd4M4J6pi9qnoauPiezXcCJ4bbJ9j8ZVm7HWqbhao6V1XPDbffAi4tMz7pvtulrrWYIuzXAr/Ycv8M81rvvYAfJHk2ycbUxWzj6JZltt4Ajk5ZzDb2XMZ7nd6zzPhs9t0iy5+P5QDd+91aVX8G/CXw10N3dZZq8xhsTnOn/wT8KZtrAJ4D/mHKYoZlxh8G/raqfrX1sSn33TZ1rWW/TRH2s8B1W+5/ZNg2C1V1dri+ADzK5mHHnJy/tILucH1h4nreVVXnq+rtqnoH+CYT7rthmfGHge9U1SPD5sn33XZ1rWu/TRH2Z4AbknwsyYeAzwCPT1DH+yS5chg4IcmVwKeY31LUjwPHh9vHgccmrOV3zGUZ752WGWfifTf58uer+mrnPb72+Q42R+T/G/jyFDXsUNefAP81XF6aujbgITa7df/H5tjGXwF/ADwFvAr8J3D1jGr7FzaX9n6BzWAdm6i2W9nsor8APD9c7ph63+1S11r2mx+XlZpwgE5qwrBLTRh2qQnDLjVh2KUmDLvUhGGXmvh/16RS7OjrAx4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(10000):\n",
    "    i, j = np.random.choice(28, 2)\n",
    "    markov_random_field[1, i, j].send_message(proprange=3)\n",
    "restored_img = np.zeros_like(noisy_img)\n",
    "for i, j in itertools.product(range(28), range(28)):\n",
    "    restored_img[i, j] = np.argmax(markov_random_field[0, i, j].proba)\n",
    "plt.imshow(restored_img, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
